Efficient algorithms for rigid body integration using optimized splitting methods and 

exact free rotational motion 
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In this note we present molecular dynamics integra- 
tion schemes that combine optimized splitting and gra- 
dient methods with exact free rotational motion for rigid 
body systems and discuss their relative merits. The al- 
gorithms analyzed here are based on symplectic, time- 
reversible schemes that conserve all relevant constants of 
the motion. It is demonstrated that although the algo- 
rithms differ in their stability due to truncation errors as- 
sociated with limited numerical precision, the optimized 
splitting methods can outperform the commonly-used ve- 
locity Verlet scheme at a level of precision typical of most 
simulations in which dynamical quantities are of interest. 
Useful guidelines for choosing the best integration scheme 
for a given level of accuracy and stability are provided. 

Hamiltonian splitting methods are an established tech- 
nique to derive stable and accurate integration schemes in 
molecular dynamicsri The strategy of these methods is to 
split the Hamiltonian of the system into parts whose evo- 
lution can be solved exactly. Using the Campbell-Baker- 
Hausdorff formula^., splitting algorithms can be presented 
as products of exactly solvable propagation steps, involv- 
ing more factors for higher-order schemes.— The resulting 
algorithms can be optimized by adjusting the form of the 
splitting to minimize error estimatesi^ 

Recently, second- and fourth-order symplectic integra- 
tion schemes for simulations of rigid body motion, based 
on the exact solution for the full kinetic (free) propagator, 
have been proposed.^ While this exact solution involves 
elliptic functions, elliptic integrals and theta functions,^ 
there exist efficient numerical routines to compute ellip- 
tic functions,^ and the computation of elliptic integrals 
and theta functions can be implemented efhciently^ or 
avoided altogether using a recursive methodi^ Employing 
the exact free rotational motion, the resulting splitting 
method leads to demonstrably more accurate dynamics 
for systems in which free motion is important;^ Further- 
more, using the exact kinetic propagator, any splitting 
scheme for integrating the dynamics of point particles 
can be transferred to rigid systems. Here we analyze the 
combination of the exact kinetic propagator and opti- 
mized splitting and gradient-lik o^'^'^° approaches. 

For a system of rigid bodies, a phase space point T 
is specified by a center of mass position q^, an attitude 
matrix Si, and translational and angular momenta pi and 
£i for each particle i of mass rrii . Given the Hamiltonian 
H = T + V, where T and V are the kinetic and potential 
energies, respectively, the time evolution of the point T in 
phase space is governed by f = {H, T} = {T, T} + {V, F}, 



in which {, } denotes the Poisson bracket. Henceforth, 
the operators {T, .} and {V,.} will be designated as A and 
B, respectively. Defining C = A + B, the solution of the 



equations of motion is formally given by T(t) 



F(0). 



While the various possible splitting schemes can be as- 
signed a theoretical efficiency,'^ the relative efficiency of 
real simulations can be somewhat different. Nonetheless, 
the estimates are useful to eliminate the least efficient 
variants. Based on our studies of second and fourth order 
methods, the most efficient integration schemes can be 
formulated using the following generic form of the split- 
ting algorithm for a single time step of size h: 



This propagator is applied t/h times to compute the time 
evolution of the system over a time interval t. Here, r] and 
^ are two real parameters, k is the order of the integration 
scheme, and e^'^ and e^'^ act on a phase space point 
r = {qt,Pi,St,li} as 

e^'T = {q, + /^p,/m„p„P,(/i)S„£J, (2) 
e^'T = {ci,,p, + hi,, S,,£, + hT,}, (3) 

where fi and Ti are the instantaneous forces and torques 
on body i, while the matrix Pi{h) propagates exactly 
Si over the time interval h in the absence of torques 
[see Ref. for specific forms for Pi(/i)]. Finally, -B(^) 
in Eq. ([T]) is a variation of B which takes the gradients of 
forces and torques into account by an advanced gradient- 
like method.'-- More precisely, the action of e^^^^'^ on a 
phase space point is given by 

e^iOhY ^ {q,,p, + h{,,S,,e, + hn}, (4) 

where the modified forces fi and torques Ti areiS 

f, = f, + Af,(e,A), T, = n + An{^,X). (5) 

The shifts in forces and torques account for commuta- 
tor corrections involving gradients. To fourth order in 
h, the shifts can be approximated by a finite difference 
approach using a small parameter A according to 



Af,(e,A) 

Ax,(e,A) 



K(q,S) 
h(q,S) 



f»(q,S)]/A, 
- x,(q,S)]/A, 



(6) 



where fi(q, S) and Xi(q, S) are the forces and torques at 
the auxiliary coordinates 

q, = q, + 2(\h%/m,, S, = Ri2i\h^ ^r^S,T,)S,. (7) 
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FIG. 1: Efficiency of integration schemes for simuiations of 
rigid water. For various values of tlie timestep h, the plot 
sliows tfie relative error versus cost (in force evaluation per 
ps) . The plots extend up to values of R where the simulations 
start to exhibit statistically-significant drift due to numerical 
round-off. The inset shows the same on a logarithmic scale. 

Here, ^i is the diagonalized moment of inertia tensor of 
the ith body [i.e. diag(/i, /2, /s)] and R(v) is the Ro- 
drigues matri»li that performs a rotation around a vec- 
tor V. Note that for ^ = 0, B{0) — B, in which case 
there are no advanced-gradient contributions. Although 
the finite difference approach introduces non-symplectic 
terms of order A^/i^, no discernible energy drift was found 
for small integration time steps h when the value of the 
parameter A was taken to be roughly 10"'*,-° 

By tuning the parameters ^ and r/, different integration 
schemes can be obtained. Choosing ^ = and 77 or 
T] = 1/2 results in the well-known second-order (fc = 2) 
Verlet scheme, in its position or velocity form, respec- 
tively. Fixing ^ = but allowing rj to vary, the prefactors 
can be minimized in front of the 0{h^) corrections, which 
gives 77 = 0.1931833275037836 as an optimal choice-^ii^ 
This scheme, which was called H0A2 in Ref. [l^, is still 
second order but is expected to be more accurate. Fi- 
nally, one can vary both 77 and ^, to make the prefac- 
tors of the 0{h?) corrections vanish to yield a fourth- 
order algorithm.— For this scheme, which we have called 
GIER4, the required values are ry = 1/6 and ^ = 1/48. 

To assess the relative computational cost of each of 
the integration schemes at a given level of accuracy, sim- 



ulations of 512 rigid water molecules using the TIP4P 
potentiali^i were carried out at liquid density of 1 g/cw? 
and a temperature of 297 K. The accuracy of the simu- 
lations was measured by calculating the ratio R of fluc- 
tuations of the total energy to the fluctuations of the po- 
tential energy at a given computational load. This load 
was estimated by using the number of force evaluations 
in a given time interval, here taken to be 1 picosecond 
(ps). At liquid densities, the computational load corre- 
lates very well with the overall CPU time since relatively 
little CPU time is required in the free motion propaga- 
tion steps. In addition, the stability of each integration 
scheme was monitored by a linear least-squared analysis 
of the drift of the total energy over a series of 10 to 50 
runs of total length 15 ps for each time step reported. 

The results of this analysis are plotted in Fig. 1, from 
which it is evident that for crude simulations requiring 
only modest energy conservation (i.e. R > 1.5%), the 
standard Verlet algorithm is the only algorithm that is 
stable. Trajectories at this level of accuracy can be used 
in sampling schemes such as hybrid Monte-Carlo. How- 
ever for R < 1.5%, arguably the upper limit of allow- 
able error in simulations from which dynamical informa- 
tion can be extracted, the optimized second-order H0A2 
scheme is roughly 1.5 times more efficient than the Verlet 
algorithm. Note that the H0A2 algorithm differs from 
the velocity Verlet scheme only in the choice of time step 
for the momenta updates and is therefore simple to im- 
plement. Interestingly, the fourth order GIER4 scheme 
is preferable if very accurate simulations are required 
{R < 0.4%) in spite of the additional computational cost 
of the modified forces and torques at auxiliary positions. 
Other fourth-order splitting schemes^ (not outlined here) 
have also been tested and found to be less efficient than 
the relatively simple GIER4. Streamlining explicit cal- 
culations of the gradients of forces and torques instead 
of utilizing finite difference methods would restore sym- 
plecticity and likely increase the value of R at which the 
GIER4 method is optimal. 
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